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SUMMARY 

The method of regularized Stokeslets (MRS) is a numerical approach using regularized fundamental 
solutions to compute the flow due to an object in a viscous fluid where inertial effects can be neglected. 
The elastic object is represented as a Lagrangian structure, exerting point forces on the fluid. The forces 
on the structure are often determined by a bending or tension model, previously calculated using finite 
difference approximations. In this paper, we study Spherical Basis Function (SBF), Radial Basis Function 
(RBF) and Lagrange-Chebyshev parametric models to represent and calculate forces on elastic structures 
that can be represented by an open curve, motivated by the study of cilia and flagella. The evaluation 
error for static open curves for the different interpolants, as well as errors for calculating normals and 
second derivatives using different types of clustered parametric nodes, are given for the case of an open 
planar curve. We determine that SBF and RBF interpolants built on clustered nodes are competitive with 
Lagrange-Chebyshev interpolants for modeling twice-differentiable open planar curves. We propose using 
SBF and RBF parametric models within the MRS for evaluating and updating the elastic structure. Results 
for open and closed elastic structures immersed in a 2D fluid are presented, showing the efficacy of the 
RBF-Stokeslets method. Copyright © 0000 John Wiley & Sons, Ltd. 
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1. INTRODUCTION 


The method of regularized Stokeslets (MRS), developed by Cortez |1|2) , is a numerical method used 
to determine fluid flow in the presence of an elastic structure at zero Reynolds number. Using the 
lineality of the Stokes Equations, the velocity field is calculated as a superposition of regularized 
fundamental solutions in which the singularity has been removed. Since microorganisms live in 
a viscosity dominated environment where the effects of inertia can be neglected, the MRS has 
been successfully used to model bacteria (e.g. E. coli and spirochetes) @0, dinoflagellates (6), 
sperm cilia ]13f|15| , biofilms |16|, filaments [17 and biflagellated algae (19). 


In each of these applications, the force that the organism exerts on the surrounding fluid is modeled 
by a regularized force term in the Stokes equation. We consider elastic structures treated as 
immersed boundaries, where the curve corresponds to a centerline approximation of a slender body 
and can be open or closed. The curves can be described in a Lagrangian form that can be discretized 
into a set of points whose trajectories wifi be calculated and updated using the MRS. These 
elastic structures generate bending or stretching forces via virtual springs or connections between 
the discretized points, which have been previously calculated by finite difference approximations. 
Alternative representations of structures via an interpolant and direct calculation of forces from 
interpolants has not been explored in the context of the MRS; however, this has been completed for 
closed curves in the Immersed Boundary (IB) method [20] 211. In this work, we are interested in 
exploring two variants of Radial Basis Function (RBF) interpolants to describe Fagrangian elastic 
structures within the MRS. 


RBF interpolation is a powerful tool for the approximation of data at scattered node locations in 
multiple dimensions. On planar domains, RBF interpolants can be viewed as a generalization of 
polynomial interpolants (which are generally guaranteed to be unisolvent only on scattered nodes in 
ID) (22) . In this work, however, we are interested in the comparison of RBF interpolants restricted 
to a circle (23] §6], (24] Ch. 17] to more standard RBF interpolants. Studies have shown that these 
restricted RBF interpolants (also known as Spherical Basis Functions or SBFs) provide spectral 
accuracy for smooth target functions (25l, and error estimates for target functions in Sobolev spaces 


are also known 126 [. To date, a complete comparison of SBF interpolants and more standard RBF 
interpolants in the parametric interpolation setting has not been completed. Of greater interest, SBF 
interpolants have been shown to be generalizations of Fourier-based interpolants (27) and exhibit 
superconvergence in the approximation of derivatives of functions on circles, spheres and tori 


In previous work |20) , we presented a parametric SBF representation of the boundaries of closed 
objects simulated in the IB method and showed that it was more accurate and less costly than the 
collection of techniques used in the traditional IB method (using piecewise quadratics to compute 
normal vectors and using finite-differences to compute forces). In addition, we showed that the 
parametric SBF model offered accuracy and convergence comparable to Fourier-based methods 
when the target shapes were infinitely smooth, and better accuracy and convergence than Fourier- 
based methods when the target shapes had only one or two underlying derivatives. In subsequent 


Copyright © 0000 John Wiley & Sons, Ltd. 
Prepared using fldauth. els 


Int. J. Numer. Meth. Fluids (0000) 
DOI: 10.1002/fld 











3 


work 0, we used the SBF geometric model within a full IB simulation of closed surfaces 
(platelets). The SBF geometric model has since been used within the Augmented Forcing method, 
a second-order forcing method for fluid-structure interaction [29j, and for forming surfaces out of 
point cloud data (on which PDEs were then solved) [30 3Tj. However, in all these scenarios, the 
SBF geometric model was used only for the modeling of closed curves and/or surfaces, and the 
basis functions themselves were always restricted to the circle/sphere. 


In this work, we explore the use of this SBF geometric model for the simulation of elastic structures, 
represented by both open and closed curves, within the MRS. Motivated by the application of 
modeling cilia and flagella, we first focus on a detailed comparison of SBF, RBF and Lagrange- 
Chebyshev geometric models for a Lagrangian structure corresponding to a planar, open curve. A 
static test case of a filament corresponding to a perturbed sinusoidal shape, a simplified centerline 
representation of a sperm flagellum, is used to investigate interpolation error as well as error in 
computing normals and second derivatives. These studies are performed on Chebyshev nodes and 
the so-called “mapped Chebyshev” nodes in parametric space, to help determine which node sets 
are well-suited for the modeling of perturbations of idealized shapes. We then present results for 
time dependent fluid-structure interaction using an SBF parametric model within the MRS. We test 
the accuracy of the SBF and RBF models within this method for a test case with an exact solution in 
comparison to finite difference approximations of forces. Results of time-dependent fluid-structure 
interaction simulations involving closed planar shapes and an open filament propagating a sinusoidal 
wave are presented. The proposed RBF-Stokeslets method compares well with the MRS using finite 
difference approximations for the forces and computation time could be potentially decreased with 
proper choice of nodes. 


The remainder of this article is organized as follows. In Section 2, we describe the adaptation of the 
SBF geometric model to the modeling of open and closed planar curves, and present a simple new 
RBF model for the same. We review the method of Regularized Stokeslets in Section 3, and describe 
how we use the SBF and RBF models within that method; specifically, we discuss modifications 
to the standard Stokeslet algorithm to incorporate the two RBF models. In Section 4, we present 
results for error on static test cases as well as time dependent simulation results for a dynamic 
elastic structure using the SBF-based regularized Stokeslets method. We conclude with a summary 
of our results, and outline possible future work in Section 5. 


2. GEOMETRIC MODELING WITH RADIAL BASIS FUNCTIONS 


In this section, we will briefly review the previously developed SBF model for closed planar curves 
and adapt the SBF representation to the modeling of open curves in 2D domains. We choose to use 
a parametric representation for these planar curves, since our target shapes (flagella, cilia, etc.) can 
naturally be represented in this fashion 


14 ]. 
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We represent a general planar curve at any time t parametrically by 

X(X,t) = (X(X,t),Y(X,t)), (1) 

where A is the parametric variable. If these curves are open and planar (and of the form Y = f(X ) 
where / is some smooth function), for a specific time t, we can simplify this parametrization to 

X(X,t) = (X,Y(X,t)). (2) 

We explicitly track a finite set of N,i points X'jit ),..., X% d (t), which we refer to as data sites. 
Here X d (t) := X(X d ,t), j = 1,..., Nd, and we refer to the parameter values Xf ,..., X% d as the 
data site nodes (or simply nodes). For a general planar curve, we construct each component of X 
by using a smooth parametric RBF interpolant of each coordinate of the data sites. 

We now explain how to construct two types of RBF interpolants to the y-coordinate Y (A, t ) using 
the data (A?, Yf(t )),..., (\% d , Y$ d (t)). Define Y(X,t) by 


N d 

Y(X,t) = ^2 c l(t)(j>{rx, k ) , (3) 

k —1 

for coefficients c k (t) and scalar-valued radial kernel <fi(r), which we discuss below. The argument 
r corresponds to the Euclidean distance between the points on the unit circle for a Spherical Basis 
Function (SBF) and a standard distance argument for an RBF interpolant. Specifically, we define r 
as follows for the SBF and RBF interpolant, 

f\,k = y^2 — 2cos(A — A£) (SBF), (4a) 

rx, k = \X-X d k \ (RBF), (4b) 

corresponding to the distance between a point A and X k . For the geometric modeling of closed and 
open planar curves, we use the multiquadric (MQ) radial kernel function, given explicitly by 

cj)(r) = s/l + (er) 2 , (5) 


where £ > 0 is called the shape parameter. Regardless of whether we use the SBF or the standard 
RBF, to have Y (A, t) interpolate the given data, we require that the coefficients c k ( t ), k = 1,..., N,/ 
for some t satisfy the following system of equations: 
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where r,-.fc will be the distance between A d and A d using either Equation j [4a| > or ( |4b[ i for the SBF or 
RBF, respectively. Since Tj.k = rkj, the matrix A is symmetric. For the MQ kernel in conjunction 
with the SBF, A is non-singular if A d is not a multiple of 27 t; the standard RBF is non-singular for 
all distinct values of A d 124 32]. Though A may be highly ill-conditioned for small values of e, the 
interpolant is well-defined even in the limit e —> 0 | 22||27| . We will compare the performance of both 
SBFs and RBFs on a static test problem and comment on selecting e in Section 4. The coefficient 
vector qf (t) can be computed using the same interpolation matrix A and setting the right hand side 
of Equation (|6]) to be the vector X d (t). We note that our choice of the multiquadric is somewhat 
arbitrary. One could use any of the infinitely-smooth RBFs in these scenarios. 


A somewhat more standard choice for this form of ID interpolation problem would be to use a 
polynomial interpolant at Chebyshev nodes. Indeed, these interpolants are known to converge at 
an exponential/geometric rate when approximating analytic (or sometimes even simply smooth) 
functions. For completion, we therefore also present a well-known stable algorithm for computing 
Fagrange interpolants at Chebyshev nodes. In this work, we will refer to these interpolants as 
Fagrange-Chebyshev interpolants, since the term “Chebyshev interpolant” could also refer to 
interpolation with a Chebyshev polynomial (which we do not consider in this work). 


Assuming a parametrization of open planar curves on an interval [a, 6], the corresponding Chebyshev 
nodes are given by 


K = 


b — ( 




cos 


2k- 1 
2 N d ' 


(7) 


for k = 1,..., N d . Note that the Chebyshev nodes do not include the endpoints of the interval. We 
must now form the polynomial interpolant at the Chebyshev nodes. It is well-known that forming 
the Vandermonde interpolation matrix corresponding to the monomials is extremely ill-conditioned 
at equispaced nodes. To ameliorate this difficulty, we form the Fagrange interpolating polynomial 
at the Chebyshev nodes. For the Y (A, t) component, this may be written as 


N d 

Y (x,t) = ^y fc d 4(A), 

k =1 


where 


4(A) 


n 

l<m<k 

m^k 




( 8 ) 


(9) 


are the Fagrange basis polynomials. To further ameliorate the ill-conditioning, we choose to use the 
barycentric form of the Fagrange interpolating polynomial. First, we define the barycentric weights 
to be 


1 

Wk -ih<><N d (K-\ d y 


( 10 ) 
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which allows us to write the Lagrange basis functions as ^(A) = 1(A) "L , where €(A) = (A — 
A d )(A — A 2 )... (A — A d - ; ). The barycentric form of the Lagrange interpolant is very efficient to use 
for a fixed parametrization, and quite stable for hundreds of points. 

As mentioned earlier, we model open planar curves by the simpler parametrization X(X,t) = 
(A, Y (A, t)). In this case, since the parametrization is a simple linear function in the x-direction, we 
only need to interpolate the y-coordinate of the data sites and compute the coefficient vector (t). 
All derivatives of the x-coordinate can be computed analytically. If one does require an interpolant 
X(A, f), we recommend replacing the scalar-valued radial kernel <f>(r) in Equation Q with the term 
| A — Af |. This is equivalent to a linear spline (which would be exact for linear functions). In our 
experiments, we also found that simply using the RBF in Equation ( |4hj i also produces reasonable 
results, with the RBF interpolant being generally better conditioned than the linear spline. 

Once we have an SBF or RBF interpolant, we may wish to evaluate the interpolant at some set of 
N s evaluation nodes, {A® }^ 1 , where N s » Nd . Further, we may wish to compute derivatives 
of the interpolant and evaluate those derivatives at either the data site nodes or the evaluation 
nodes. We will now describe how to perform these operations without explicitly computing the 
coefficient vectors c d (t) and c d (t). Again, for simplicity, we focus on evaluating and differentiating 
the interpolant Y(A, t). Formally, from Equation (J6]», we have 


c Y (t)=A~ 1 Y d (t). (11) 

We call the Cartesian points corresponding to the evaluation nodes the sample sites. Now. let 
Y s (t) be the vector with N s entries that contains the y-coordinates of the sample sites, obtained 
by evaluating the interpolant Y( A, t) at the set of evaluation nodes. To obtain the vector Y s (t), we 
perform the following operation: 


Y s (t) = Bc Y d {t) = BA~ 1 Y d {t), ( 12 ) 

where {for j = 1,2,..., N s , k = 1,2,..., N d , and r J:fe is defined by Equation ( |4li| ) 
or ( |4hj i for the SBF and RBF, respectively. E = BA~ l is the time-independent N s x N r j evaluation 
matrix. One can obtain the x-coordinates of the sample sites similarly by the operation X s (t) = 
EX d (t). Defining the matrices X d (t) = [X d (t ) Y d (f)] and x\t) = [X s (t) Y s (t)\, then J S (f) = 
EX (t). For open planai' curves, we find that it is sufficiently accurate to use a single evaluation 
matrix E, though one could also take the approach of defining two evaluation matrices (one 
corresponding to a piecewise-linear interpolant for the x-coordinate and the other corresponding 
to an SBF or RBF interpolant for the y-coordinate). 

In order to compute quantities like tangents, normals and forces, we also need to be able to compute 
derivatives of the coordinate interpolants with respect to the parameter A and evaluate them at either 
the data site nodes A d or the evaluation nodes A s . To compute the n th - order derivative of either an 
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SBF or RBF interpolant and evaluate the derivative at the data site nodes, we define the matrix 




, j,k = l,...,N d , 


A=A? 


( 13 ) 


which is an N d x N d matrix. To compute the -order derivative of an RBF interpolant and 
evaluate the derivative at the evaluation nodes, we similarly define the matrix where 

j = 1,..., N s , and k = 1,..., N d . If we explicitly computed the coefficient vectors c d (t) and ( t ), 
the products of the matrices B'f, and with those vectors would yield derivatives of interpolants 
at the data sites nodes and evaluation nodes respectively. Flowever, since it is cheaper to use a 
coefficient-free formulation, we will now define analogs to Equation ( fl2j ) that perform coefficient- 
free differentiation. We thus define the following differentiation matrices: 






(14) 


which yield derivatives of the A’(A. t) and Y(X,t) interpolants at data site nodes and sample site 
nodes when post-multiplied by the vectors A d (i) and Y d (t). These differentiation matrices thus 
combine the operations of differentiation and evaluation of the interpolants. 


It is important to note that since the evaluation and differentiation matrices implicitly contain the 
interpolation operation (because they contain the term A -1 ), these matrices can be applied to any 
-/V^-long vector of function samples at the data sites to obtain either samples of the function at 
sample sites, or samples of derivatives of the function at data sites or sample sites. For a more 
detailed description of these interpolation and differentation matrices, we refer the reader to previous 
work (21]. 

The barycentric form allows for fast evaluation of the Lagrange-Chebyshev interpolant. It also 
allows one to compute derivatives quite easily at the data site nodes. Flowever, it is not clear how 
to evaluate the derivatives at the sample site nodes using this form of the interpolant, unlike with 
RBF/SBF or other interpolants. Flowever, we use a simple alternative to evaluate derivatives. We 
first use the barycentric form of the interpolant to obtain derivatives of the interpolant at the data site 
nodes; then, we once again use the barycentric form of the interpolant to interpolate the derivative 
samples themselves, and evaluate that interpolant at the sample site nodes. For polynomials, we 
believe this process is quite reasonable, with the slight additional algorithmic difficulty being more 
than compensated by the stability, efficiency and flexibility of the Lagrange-Chebyshev interpolant 
in barycentric form. With this framework in place, we will compare the Lagrange-Chebyshev 
interpolants to the SBLs and RBLs in Section 4.1. 

Lor all our tests on open planar curves, we use nodes contained in [0,1], simply to facilitate 
the kinds of parametrizations used within the MRS. It is important to find a set of interpolation 
nodes Xf ,..., A d - rf that results in a stable interpolant. When the shape of interest is a closed curve 
homeomorphic to the circle, evenly-spaced values indeed give us stable interpolants. When moving 
to an open curve, for stability, we must now choose between different types of clustered nodes 
mapped to the interval of interest. Regardless of the choice of interpolation nodes, we choose to 
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always evaluate our interpolant at a set of equispaced nodes in the interval for static tests. In full 
Stokeslet simulations involving open curves, we track and interpolate the data sites using a fixed 
parametric interpolation matrix built out of clustered nodes. The reasoning for this will become 
apparent when the results of the static tests for different types of nodes are seen. We will explore 
the effects of different node choices in Section ITT! 


3. USING THE RBF MODEL WITHIN THE REGULARIZED STOKESLET METHOD 


In this section, we will first review the method of Regularized Stokeslets. We will then discuss how 
both the RBF representations are used (in lieu of finite difference approximations) for geometric 
modeling within this method. 

Consider the scenario where an elastic curve (either open or closed) is immersed in a viscous, 
incompressible fluid described by the Stokes equations. In two or three dimensions, these equations 
are as follows: 


HA u = Vp-f, (15) 

V • u = 0, (16) 

where /y is the viscosity of the fluid, p is the pressure, u is the velocity and / is a body force acting 
on the fluid. For a single point force, the exact solution for the resulting flow is a fundamental 
solution called a Stokeslet. The Stokeslet can be used to evaluate the fluid velocity at any point off 
of the structure and is singular if evaluated at the location of the point force. Thus, Cortez derived 
the method of regularized Stokeslets, removing the singularity by regularizing the forces |1]|2j. This 
method can also handle multiple point forces since the Stokes equations are linear and the solution 
can be obtained by a superposition of regularized fundamental solutions. 

The elastic curve that exerts forces on the surrounding fluid is typically represented by a Lagrangian 
description, X(A, t) with parameter A at time t. In the MRS, the continuous formulation of a body 
force / on a closed or open curve T is written as 

/(*)= J F{X,t)(f>s(x - X)d\, (17) 

where F is the force per unit length generated along the structure, x is any point in the fluid domain 
(including possibly a point on the structure), and <j >,5 is a radially symmetric smooth approximation 
to a delta function. In this formulation, the fluid feels the structure primarily in a region around 
the curve, governed by the regularization parameter S. We note that in practice, the Lagrangian 
structure is discretized into a set of points and the body force / is found by approximating the 
integral in Equation fT7} with a quadrature rule. 
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Using a regularized body force as in Equation the regularized solutions for the pressure 
and velocity at the IB points (and indeed, in all of space) can be derived. If there are N s forces 
T k = —F k A A at IB points X k that are equally spaced by AA, the pressure and velocity in R 2 are 
explicitly given by 


n 3 

p(x) = Y / (F k -VG 5 (x-X k )), (18) 

k =1 

1 N * \ 
u{x) = - Y, ( + (^fc • V)VB { (x - X k ) - T k G s (x - X k ) J . (19) 

" k —i ' ^ ' 


Gf, is a smooth approximation of the Green’s function and Bs is a smooth solution to the 
biharmonic equation such that AGs = <t>s and A Bs = Gs- Both Gs and Bs are assumed to be 
radially symmetric, simplifying expressions such as XBg(x — X) = B' s (r)(x — X)/r where r is 
the Euclidean distance between x and X. In R 2 , the blob function cj>s{x — X) = 2 Tr^+s 2 ) 5 / 2 ^ as 
corresponding regularized functions Gs and /j': 


Gs(r) 


B's(r) 


In {R + S)-^ , 

2 

2r In (R + 6) — r — — -- 

It + o 


( 20 ) 

( 21 ) 


where R = y/r 2 + 5 2 . We remark that taking the limit 5 —> 0 for all these expressions will yield 
the traditional singular expressions for velocity and pressure in the standard Stokeslet formulation. 
Equations (fTRI (1~9|) are an exact solution for the given regularized forces in Equation ([17} and the 
velocity field in Equation (17} is everywhere incompressible. For more detail, we refer the reader 
to 0. 


In the application of ciliary or flagellar movement, an elastic structure is moving and generating 
forces in a time dependent manner. At each time step, we will have a new body force f in the Stokes 
equations, (T5}-©. As a result of the new force, we can determine the new pressure and velocity 
field using the MRS and evaluating Equation (T8}-(l9}. In the MRS, similar to the IB method [ f33~| , 
we will assume that the immersed elastic structure is moving at the local fluid velocity. In a standard 
MRS algorithm, one would evaluate Equations (T8}-(T9} at each of the N s IB points and then update 
their location by solving a system of 2 N s differential equations of the form u {X) = . 

When using either RBF representation within the method of regularized Stokeslets, the first step 
is to use the data sites, the appropriate force model and the operators from the previous section to 
generate forces at the sample sites. We then compute the pressures and velocities at the N d data sites 
using the forces at the N s sample sites. In other words, we use Equations (18} and (19} but evaluate 
those equations only at the locations x = the set of data sites. While the number of terms 

(TV, 5 ) in the summation for the velocity remains the same, we now only need to compute N d such 
summations, where N d « N s . By using N s terms in each of the N d summations, we maintain the 
accuracy of the quadrature. The RBF representations thus offer a slight reduction in computational 
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cost, but a great improvement in accuracy for a given computational cost when compared to the 
traditional algorithm for the method of regularized Stokeslets. 

We outline the RBF-Stokeslets algorithm below. 


Given a set of data sites X (t), we use either RBF interpolant to generate a set of sample sites 
X (t) (at equispaced parametric nodes) using a variation of Equation ( fl2) . 

Evaluate forces at the sample sites X (t) using the differentiation matrices of the RBF 
interpolant as outlined in Equation {74]). This will vary based on whether forces depend on 
tangents, normals, etc., and may involve intermediate steps where derivatives of quantities are 
computed at data sites. 

Evaluate the velocity in Equation {l9| at N d data sites. That is, evaluate ufX'j) for j = 
1,..., N d where the summation is over the N s sample sites where N d <C N s . 

Update the location of the N d data sites by solving a system of 2 N d differential equations, 
u (Xj) = for j = 1,..., N d . 


After this series of steps, the same process repeats at the next time step. We note again the 
importance of two sets of sites, data and sampling for the SBF/RBF version of the MRS. The 
SBF/RBF parametric model of the structure allows efficient and accurate force calculations at N s 
sample sites (evenly spaced in the parametric variable A). Additionally, the computational expense 
of summing over point forces in Equation {T9) l to update the location of the IB points is decreased 
by only calculating velocities at N d data sites. In the results presented in Section 4 for the dynamic 
elastic filaments, we will use a simple Euler method for the calculation of the updated location of 
the IB data points. 


4. RESULTS 


4.1. Results for static open planar curves 


In (20), we tested an SBF interpolant against Fourier-based interpolants on two types of closed 
objects: a shape represented by an infinitely-smooth function and another shape represented by a 
twice-differentiable function. In those tests, the SBFs matched the accuracy and convergence of the 
Fourier interpolants on the smooth object, and were more accurate than Fourier interpolants for the 
object represented by the twice-differentiable function. 

Our goal in this section is to compare the ability of the SBF to approximate twice-differentiable 
functions against the ability of the standard RBF interpolant to do the same. In addition, we study 
the behavior of both interpolants against the Lagrange-Chebyshev interpolant. It is important to 
note that both the RBF and SBF interpolants generalize trivially to 2D parametric interpolation on 
scattered nodes, while such a generalization of the Lagrange-Chebyshev interpolant is non-trivial. 
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We start with a simple test function of the form 


X id eai{ A) = ( X idea i(\),Y idea i(\ )) = (A, 6sin(27rA — w)), A G [0,1]. 


( 22 ) 


This test function represents an idealized sinusoidal shape, and is often the parametrization used for 
a sperm flagellum within regularized Stokeslet simulations |7]|8 TUJ 1~2| . Obviously, if this shape 
always stayed sinusoidal, it would be impossible to have a better interpolant than a parametric 
Fourier interpolant, which contains a sine term in its expansion and can consequently reproduce the 
sinusoidal wave with spectral accuracy. However, over the course of a simulation, it is likely that 
the immersed elastic curve may deviate from the ideal sinusoidal shape. The resulting shape may 
be represented either by a smooth function or by a finitely-differentiable function, depending on the 
nature of the spatial discretization of the regularized Stokeslet method. We now consider the case of 
twice-differentiable functions. 


Roughly-perturbed Curve 



Figure 1. Test object used for error analysis. The ideal shape is shown with a solid line (-) and the roughly 
perturbed version is shown with a dashed line (- -). Note the visual similarity of the perturbed shape to the 
ideal shape despite the difference in their mathematical representations. 


We define a perturbation of this idealized shape (see Figure [TJ. This perturbation results in a curve 
that, while visually not very different from the idealized shape, has only two underlying derivatives. 
We only perturb the y-components of the shape for this test. The shape is given by: 


Y P 


1.0 + zlexp 


(l — cos 2 (27tA)) 1 ’ 5 '^ 


Yideal • 


(23) 


In the above equation, A controls the amplitude of the perturbation, with a controlling the width of 
the perturbation. For the tests that follow, we use A = 0.04 and a = 0.9. This allows us to reproduce 
the conditions of the test in [20], now allowing us to compare alternatives to Fourier interpolants. 
For all our tests, we set b = 0.05 in Equation ( |22j ). These are the values used for the plot in Figure [T] 


Since we are parametrically interpolating open curves, it is unwise to use equispaced nodes lest we 
encounter the Runge phenomenon. We therefore present results only for clustered nodes. However, 


Copyright © 0000 John Wiley & Sons, Ltd. 
Prepared using fldauth. els 


Int. J. Numer. Meth. Fluids (0000) 
DOI: 10.1002/fld 












12 


since Fourier interpolants are likely unstable on clustered nodes, and are inappropriate for the 
modeling of curves represented by twice-differentiable functions, we do not present any results 
for this case. In our experiments, the poor accuracy of Fourier interpolants in interpolating such 


shapes was quite similar to that seen in [201. 


Before we proceed, we quickly remark on the choice of interpolant for infinitely-smooth open 
curves. For curves that are represented by periodic and smooth functions, the obvious victor is 
the Fourier interpolant, just as in the closed curve case seen in pO) . However, similar to the results 
from that work, SBF interpolants are competitive with the Fourier interpolants in this scenario, with 
the only drawback being that errors cannot be pushed all the way down to machine precision for 
the SBF interpolants due to ill-conditioning. This may be somewhat surprising, since the use of 
SBFs implies that we are parametrizing open curves on the circle. However, we have found that 
this approach works perfectly well. Of course, since the curves are open, SBFs have no advantage 
over standard RBFs (unlike in the closed case, where RBFs fail to capture the periodic point well). 
However, SBFs appeal' to be equivalent to RBFs in this setting. Depending on the performance of 
SBF interpolants in the following sections, this implies that one could write very general software 
based on SBF interpolants that would work excellently for the modeling of both closed and open 
curves. We will not present results for the recovery of infinitely-smooth functions corresponding 
to open curves with either SBFs or RBFs in this work, and focus on the more interesting case of 
twice-differentiable functions representing open curves. 


For the modeling of open, twice-differentiable curves with RBF and SBF interpolants, we use 
two types of clustered parametric nodes: standard Chebyshev nodes, and the so-called “mapped 
Chebyshev nodes” introduced by Kosloff and Tal-Ezer. For brevity, we will refer to the mapped 
Chebyshev nodes as KTE nodes. In |34j, it was shown that the Lebesque constants corresponding 
to RBF interpolation at the KTE nodes grow much more slowly than in the case of RBF 
interpolation at Chebyshev nodes. For a parametric Chebyshev node \ cheb defined in Equation ([7|, 
the corresponding KTE node is given by 


sin _1 (aA cftei ') 
sin -1 (a) 


(24) 


When a = 1, the nodes are equispaced, while the limit a —> 0 recovers the Chebyshev nodes. 


4.1.1. Interpolation at Chebyshev nodes We first compare the ability of the SBF, RBF and 
barycentric Lagrange interpolants to approximate Yp, its normals and its second derivatives with 
respect to A. All interpolants are built on the same set of Chebyshev data sites nodes in (0,1), and 
all interpolants are evaluated at a set of N s = 400 equispaced sample site nodes in [0,1]. The results 
are shown in Figure [2] For each case in Section 4.1, the error is defined as the maximum of the 
point-wise C 2 errors for each of the sample sites. 

From Figure [2j we see that both the SBF and the RBF have very similar errors in all cases for 
the same shape parameter at the Chebyshev nodes. The figure on the top left shows that both 
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Comparison of Evaluation Errors 


O SBF-Chebyshev Interpolant (e = 7.0) 
* 


Comparison of Errors in Normals 


20 40 60 

Number of data sites 


20 40 60 

Number of data sites 


O 

* 

o 


SBF-Chebyshev Interpolant (£ = 7,0) 
RBF-Chebyshev Interpolant (£ = 7.0) 
Lagrange-Chebyshev Interpolant 



Figure 2. Comparison of interpolants at Chebyshev nodes when interpolating and evaluating a shape 
represented by a twice-differentiable function (top left), its unit normals (top right) and its second derivative 

(bottom). The shape is given by Equation l|23|l. 


the SBF and RBF interpolants have errors that are quite similar to the errors of the barycentric 
Lagrange interpolant, albeit slightly lower for higher values of N d . However, the figure on the top 
right and the one on the bottom show that the SBF and RBF interpolants have lower errors when 
computing normals, and much lower errors for the second derivatives. This may be a consequence 
of the difficulty in interpolating a function and evaluating derivatives at a different set of nodes in 
the barycentric formulation of the Lagrange-Chebyshev interpolant. Interestingly, the SBF seems to 
have very similar errors to the standard RBF, despite the difference in the distance argument passed 
in to the basis function </>(r). While SBFs are preferred for the interpolation of analytic and/or very 
smooth periodic functions corresponding to closed curves (as seen in (20)), this indicates that they 
can safely be used in the context of modeling open, twice-differentiable shapes with high accuracy 
as well. 
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4.1.2. Interpolation at KTE nodes We now compare the SBF and RBF interpolants built on KTE 
nodes against the barycentric Lagrange-Chebyshev interpolant in the approximation of Yp, its 
normals and its second derivatives. The KTE nodes were generated with a = 0.85 and a shape 
parameter of e = 7.0. Note that at this value of a, one cannot safely use the Lagrange interpolant 
with the KTE nodes, since these nodes are quite close to evenly-spaced. The results are shown in 
Figure [3] 




O SBF-KTE Interpolant (c = 7.0) 

* RBF-KTE Interpolant (c = 7.0) 

C> Lagrange-Chebyshev Interpolant 


Comparison of Errors in Normals 


20 40 60 

Number of data sites 



Figure 3. Comparison of SBF/RBF interpolants at KTE nodes, with the Lagrange interpolant at Chebyshev 
nodes, when interpolating and evaluating a shape represented by a twice-differentiable function (top left), 
its unit normals (top right) and its second derivative (bottom). The shape is given by Equation (|23)l. 


A careful examination of Figure [3] shows that both the SBF and RBF interpolants have even 
lower errors at KTE nodes than they do at Chebyshev nodes. Consequently, they have much lower 
eiTors than the barycentric Lagrange-Chebyshev interpolant for evaluation, normal computation and 
computation of second derivatives. Once again, the SBF and RBF interpolants have almost identical 
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errors for a given value of N d when interpolating this twice-differentiable function. This shows the 
promise of the SBF in interpolating functions of very different smoothness. 


We note that even for infinitely-smooth target functions (results not shown), the SBF-KTE 
interpolant can beat both the SBF interpolant at equispaced nodes and the RBF interpolant at either 
set of nodes. For some finite and small value of Nd, the SBF-KTE interpolant actually has lower 
errors than the Fourier interpolant (the latter collocated at an equispaced set of nodes). The errors 
for the Fourier interpolant, unlike those for the SBF-KTE interpolant, can be pushed all the way 
down to machine precision. To do so for the SBF interpolants, we would require a so-called “flat” 


algorithm (these tend to be more expensive), such as the RBF-QR algorithm 1271. 


Interestingly, in our experiments with the RBF-QR algorithm, we have observed that as e —> 0, the 
SBF-KTE interpolant becomes unstable, and a change to the standard equispaced nodes is required; 
similarly, the RBF-KTE interpolant must be swapped out for the RBF-Chebyshev interpolant as 
e —> 0. Both these changes are to be expected, since in that limit, the SBF and RBF interpolants 
recover the Fourier and Fagrange interpolants respectively. We do not explore such strategies here, 
as our goal is to attain a cost-effective and highly accurate approximation of both function and 
derivative approximation for small values of Nd and non-zero values of s. 


4.1.3. Interpolation at KTE nodes with “best” shape parameters Until now, we have made the 
rather simple choice of selecting a single shape parameter for all values of N d for both the SBF and 
RBF interpolants. This choice will negatively affect the errors for smaller values of Nd- Indeed, for 
a given value of N d , the correct strategy is to pick a single shape parameter that results in the lowest 
error. It is possible to do so for each value of Nd in an error plot. We run such an experiment at the 
KTE nodes (which we have already determined to be a better node set than the Chebyshev nodes 
for the SBF and RBF) using a = 0.85. The results of this test are shown in Figure |4] 

As expected, both the SBF and RBF interpolants now have lower errors for the smaller values of N,i 
than they did when using the fixed value of e = 7.0. Furthermore, while the error curves do not look 
quite as pleasing as those in Figure[3] it is clear that even with these optimized shape parameters, the 
SBF and RBF interpolants converge reasonably well. In addition, the errors for the SBF and RBF 
interpolants are now lower than those of the Fagrange-Chebyshev interpolant even for small values 
of Nd, which was not the case for the studies with fixed e. While there exist a few sophisticated 
algorithms to optimize the shape parameter, we simply pick the best ones by sampling densely in 
e-space, and computing errors for all those values of e for a given number of data sites N d . For 
reproducibility, we document the results in Table [I] 

While the values in Table [T] look highly tuned, this is simply an artifact of the range of samples 
used to generate the different values of e. For a given Nd, one can certainly use values of the shape 
parameter that are close to but different from those presented in the table, with little change in 
the error. This table is merely intended as a guide to selecting the shape parameter in preparation 
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Figure 4. Errors in evaluation, normal computation and calculation of second derivatives for the SBF (left) 
and RBF (right) interpolants for carefully selected shape parameters. 


N d 

£sbf 

£rbf 

8 

2.5 

2.6 

16 

3.2 

3.2 

24 

3.0 

2.9 

32 

3.8 

3.9 

40 

3.8 

3.6 

48 

4.5 

4.7 

56 

5.7 

5.9 

64 

8.2 

8.0 

72 

8.9 

OC 

OC 

80 

9.0 

8.6 


Table I. “Best” values of the shape parameter e corresponding to different numbers of data sites N d for the 
SBF and RBF interpolants. 100 samples of e were used as the candidate set. 


for a full time-dependent simulation, and for demonstrating the value in selecting smaller shape 
parameters for smaller N d . 


4.2. Stokeslet simulations on closed shapes 


In this section, we show results for simulations of closed elastic objects immersed in a viscous, 
initially stationary fluid. Error analysis of RBFs with closed curves has been previously detailed 
in the context of the IB method |20| , thus our results focus on showing that this method works 
for closed curves in the context of the regularized Stokeslet method. First, we present results for 
a stationary problem that has an exact solution. We then show results for a time-dependent fluid- 
structure interaction simulation. For both test cases, we use the multiquadric (MQ) radial kernel 
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function given in Equation Q and the SBF given in Equation ( |4a| ) to interpolate the y-coordinate at 
the data sites. 


4.2.1. Computation of tangential forces This test problem was presented in [1] (example 46). The 
forces F on the an elastic structure corresponding to the unit circle are prescribed to be purely 
tangential to the boundary, 

F(A) = -2sin(3A)^^ (25) 

for parametric variable A. The unit circle and corresponding tangential forces (black arrows on 
the structure) are shown on the left plot in Figure [5] The exact solution for Stokes equation with 
viscosity p, = 1, as detailed in [TJ and shown on the right in Figure [5] has a discontinuous pressure 
and velocity gradient (Vp and Vw) across the boundary of the elastic structure and continuous p 
and u everywhere in the infinite 2D fluid. 



Figure 5. Exact solution for a unit circle with tangential forces, presented in |TJ. The left figure shows the 
structure, the unit circle, with forces on the structure shown with black arrows. Iflie line of arrows correspond 
to evaluating the fluid velocity at the points on the line (a;, 1/5). On the right, the exact solution for the pressure 
p and velocity u = (u, v) evaluated at y = 1/5 for the x values displayed. 

We use Nd = 25 equispaced data sites for the representation of the circle with A € [0, 27 t). The 
derivative matrix in Equation l p~4} is used to compute the tangents <9X(A )/OX on the circle in order 
to calculate the force in Equation ( [25] ). We evaluate the forces at N s = 400 equispaced sample sites, 
and then compute the pressure and velocity using Equations (fTTsji ([To]) at a set of markers lying 
on the line (x, l). 0.1 <x< 1.8. We compare the errors in the solution obtained from the SBF 
representation with those in the solution obtained by using the analytically correct tangents on the 
circle. The error reported is the absolute value of the difference of the solutions obtained by using 
the two methods. The results are shown in Figure [ 6 ] 

The results in Figure [ 6 ] show that despite using only N,/ = 25 data sites to represent the structure 
with SBFs and sampling the forces at N s = 400 sample sites, the absolute value of the error in the 
SBF-based Stokeslet approximations to the pressure and velocity fields are almost identical to those 
obtained when using analytically computed tangents at N s = 400 IB points using MRS. For both 
methods, the error decreases as you move away from the boundary at x = 1. Similar results were 
also observed for the RBF-based Stokeslet approximations. On the circle, one can expect similar 
performance from a Fourier interpolant, though the Fourier interpolant will offer lower accuracy on 
roughly-perturbed closed shapes for the same number of data sites | 20 ) . 
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Error in Pressure 




Figure 6. Errors in the pressure and velocity fields when using the SBF model. The figure on the left shows 
the errors in the pressure for the test case presented in |li|, while the figure on the right show the errors 
corresponding to each component of the velocity field. A Q = 25 data sites and N s = 400 sample sites are 
used for the SBF representation and 400 sites are used for the exact force evaluation with the MRS. The 
numerical solutions are compared against the exact solution given in |Tj. 



Figure 7. Errors in the pressure and velocity fields when using finite difference approximations to the 
tangents. The figure on the left shows the errors in the pressure for the test case presented in QJ, while 
the figure on the right show the errors corresponding to each component of the velocity field. N s = 800 IB 
points are used for the finite differences. The numerical solutions are compared against the exact solution 

given in |Tj . 


To understand the power of global representations, we repeat the test with a second-order finite 
difference approximation to the tangent. In the Stokeslet and IB literature, finite difference 
approximations to derivatives at IB points are often computed using differences between points 
on the boundary that are halfway between the IB points ( e.g (35)). This is equivalent to using 
N s = 800 IB points (or N s = 400 IB points with “tight” finite difference stencils). Figure [7] 
shows the results for the pressure and velocity obtained with the finite difference approximation 
at N s = 800 IB points. Clearly, the error for the pressure is significantly higher when using the 
finite difference approximations to the tangents, in contrast to the results achieved using N,i = 25 
data sites and N s = 400 sample sites with the SBF representation. On the right, the finite difference 
approximations at N s = 800 leads to fairly low error, but it oscillates dramatically and has mean 
square error that is greater than using the exact formula for the tangents. This highlights a hazard 
of using finite differences to approximate geometric quantities. In comparison, the SBF solution is 
smooth and oscillation-free. Additionally, we emphasize that we are able to achieve a lower error in 
the SBF-stokeslets method than the FD method when using only Nd = 25 data sites and N s = 400 
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sample sites. We remark that (11 employs correction terms to the method of regularized Stokeslets 
on this test case to improve its accuracy. We do not use that approach in this current work. 


4.2.2. Time-dependent simulations. We will use a closed curve as our first test case, as previously 
used in other methods ]3hl|37| . We will start with a closed, elastic structure that is initialized as a 
perturbation of a circle and immersed in the fluid. The initial configuration X(A, 0) is defined as 

X(X, 0) = ((1 + /?cos(i/A)) cos(A), (1 + /3 cos(i'A)) sin(A)) (26) 

for A € [0, 27 t) with v and f3 as parameters. For this closed case, we will use Nd = 25 data sites and 
N s = 50 sample sites; each of these will be equally spaced points in the parameter A. The force will 
be proportional to the curvature at a given time point and defined as follows, 

F(A,f) = - ^«(A,t)n(A,i) (c{t) - ^ , (27) 

where k(A, t) is a signed curvature, h is the unit normal, and C[t) is the arclength of the structure at 
time t. These forces will tend to restore the shape of the closed elastic structure to that of a circle. 


t=0 t=i 



Figure 8. Results for a closed elastic structure initialized with /3 = 0.3, u = 3, Nd = 25, N s = 50, At = 
10 -3 , S = 4n/N s , and e = 1.1. Upper left is the initial configuration and the other three panels show the 
flow field with the black arrows and the color bar corresponds to the pressure. 

To compute the time dependent movement of the elastic structure, we will use Equation to 
evaluate the derivatives of the SBF y-coordinate interpolant at the data sites. This will allow for the 
evaluation of the curvature and normal vector. In the MRS, a no-slip condition 

,w 

u(x) = (28) 
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is used to update the location of the structure. Thus, as detailed in Section 3, we need to know 
the velocity at some number of points along the structure. Instead of evaluating the velocity at all 
N s = 50 sample sites, we evaluate the velocity in Equation ( fT9] ) at N d points summing over the N s 
sample sites. The location of the structure is then updated at the N d points X (t) using the the no 
slip condition given above in Equation ( |28) . This completes a single time step. 

Representative results for a closed elastic structure immersed in an incompressible fluid with 
viscosity p, = 1 are shown in Figure [8] In the upper left panel, this is the initial configuration of the 
closed circle. We show results at time t = 1,7, and 10 in the other three panels. Using Equation ( | 1 9[ i. 
we can also solve for the flow field anywhere in the 2D fluid domain. The arrows in Figure [8] 
correspond to the flow field. The background color and colorbar corresponds to the pressure as 
evaluated in Equation ( fT8j ). Although not shown, we note that the RBF version of the MRS and the 
standard MRS (forces calculated using a finite difference approximation) both compare well to the 
SBF version of the MRS results as shown in Figure [8] 


4.3. Stokeslet simulations on open shapes 


4.3.1. Computation of Tangential forces Similar to Section 4.2.1 we study the same example with 
forces F purely tangential to the elastic structure as defined in Equation ( |25T >. Here, the open curve 
is given by X(A) = (A, sin(A)) for A such that 0 < A < 27 t and we use Chebyshev and KTE nodes 
for the data sites. For this test case, we do not have an analytical solution for all points in the 2D 
fluid. However, the fundamental solution to the Stokes equation for point forces (not regularized) 
is given by the Stokeslet, which is defined at all points in the fluid not on the structure. We will 
use this as our exact solution at a set of marker points off of the structure to compare the error 
using different nodes. In Figure [9] the structure and tangential forces along the structure are shown. 
The fluid velocity is also determined along the line y = 1/5 with viscosity set to p = 1. The exact 
solution for the pressure and velocity along this line are shown in the right of Figure [9] 



Figure 9. Open curve on the left with tangential forces along the structure shown. The velocity at fluid 
markers along the line y = 1/5 are also shown. On the right, exact solution for the fluid velocity at the 

points y = 1/5 for the x values displayed. 
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Error in Velocity 

-u (RBF tangents - Chebyshev) 

- v (RBF tangents - Chebyshev) 

* u (RBF tangents - KTE) 


* v (RBF tangents - KTE) 



0.5 1 1.5 2 2.5 


Figure 10. Error of velocity at a line of marker points at y = 1/5 for the open curve with tangential forces. 
The error for the pressure and velocity are shown using shape parameter e = 1.1 for both Chebyshev and 

KTE data nodes with a = 0.85. 


We next wanted to compare the error when evaluating the velocity and pressure when using the exact 
tangential forces and when using the RBF and SBF with the multiquadric kernel (for interpolation of 
the y-coordinate) to calculate tangential forces on the structure. In Figure [TO] plots of the absolute 
difference between the solutions at each of the x-coordinate markers are shown for pressure, u 
component of the velocity, and v component of the velocity. We use N d = 50 data sites (Chebyshev 
or KTE) to compute forces and evaluate the tangential forces at N s = 200 equispaced evaluation 
nodes. The velocity and pressure are then calculated using the fundamental solution with forces at 
the evaluation nodes. The absolute error between the fundamental solution with exact tangents and 
RBF and SBF regularized Stokeslets solutions for the pressure and velocity are shown in Figure 


10 Using a shape parameter e = 1.1 for both the SBF and RBF, the KTE nodes (a = 0.85) has 


a decreased error in the approximation of the pressure and both components of the velocity in 
comparison to the Chebyshev nodes. For this a and e chosen, the SBF has less error than the RBF at 


the maker points evaluated. Similar to the case of the closed curve in Section 4.2.1 the error in the 
pressure and velocity when using finite difference methods to calculate the force is several orders of 
magnitude larger when using N d = 800 points (results not shown). 


4.3.2. Time-dependent simulations We will use a test case for an open curve corresponding to a 
filament propagating a planar sinusoidal wave of bending in a 2D fluid. Filaments propagating 
planar undulations have been studied as a simple model for a sperm flagellum |[8][38j. The open 
filament will be initialized as X(X k ,0) = (X k , 6sin(27rA)) where X k for k = 1,... ,N S are KTE 
nodes on (0,1), b is the amplitude of the sine wave, and the period of the sine wave is 1. In this test 
case, we want the open filament to bend or have a time dependent curvature corresponding to the 
sine wave 6sin(27rA — uit) with angular frequency w. The force on the structure will be determined 
as a sum of a tension force F 1 and a bending force F B defined as, 

F T =^- x {S t (\\t\\-1)t), 

B= (d^x dx^\ 

Sb (v dx 4 5 A 4 ) ’ 


(29) 

(30) 


Copyright © 0000 John Wiley & Sons, Ltd. 
Prepared using fldauth. els 


Int. J. Numer. Meth. Fluids (0000) 
DOI: 10.1002/fld 





























22 


where r is the tangent vector and f is the unit tangent. In Equation ( [29] ), this is a tension 
force that acts like an inextensibility constraint. The coefficients St and Sb are the tensile and 
bending stiffness for the elastic filament, respectively. The ideal or target shape is a propagating 
wave, X ; = (A, 6 sin(fcA — wi)). We note that both of these forces, common to immersed boundary 
methods, have already been tested and implemented for a closed curve in the RBF-Immersed 
Boundary framework 121 ]. Similar to the above description for the closed SBF - stokeslets method 
in Section |4.2.2| for a closed curve, we can calculate these derivatives using differentiation matrices 
in Equation ( |T4j i to determine the forces at each of the N s sample sites. We note that we will 
use Nd = 20 KTE data sites with a = 0.85 and N s = 40 equally spaced sample sites. For this test 
case, we use the multiquadric radial kernel in Equation ((5]) and SBF given in Equation ( |4aj i for 
interpolation of the y-coordinate. 



Figure 11. Results for an open filament initialized with Nd = 20, N s = 40, At = 5 x 10 -4 , <5 = 2/N s , 
St = 0.001, Sb = 0.1, e = 1.5, and a = 0.85 for KTE nodes. On the left, we set b = 0.01 and ui = —2n. 

On the right, we set b = 0.005 and u> = —47r. 


Results showing the location of the open filament at different points are shown in Figure [IT] We 
solve the Stokes equations using viscosity /r = 1 and the no-slip condition is used to update the 
location of the filament. Since the bending force in Equation ( |3Q| ) has a time dependent ideal or 
target location, the difference between the fourth derivative of the target location and the actual 
location of the structure will cause the structure to bend. On the left side of Figure [TT| we set the 
amplitude of the ideal structure to be b = 0.01 and the frequency to be w = — 27t. As shown from 
the time course on the left, the filament is moving and interacting with the fluid. The forces are 
driving the structure to achieve an amplitude that is approximately 0.01 and it propagates a wave in 
one unit of time, as expected. On the right side of Figure [TT[ we have decreased the amplitude of 
the ideal structure to be b = 0.005 and increased uj to -47r. As can be seen, the filament is moving 
with a smaller amplitude and is bending at a higher rate, corresponding to propagating a wave in 0.5 
units of time. In Figure [IT] the structure is plotted at the KTE nodes. We note that the KTE results 
compare well with results obtained using Chebyshev nodes with the SBF as well as using KTE or 
Chebyshev nodes with the RBF. 
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5. SUMMARY 


The method of Regularized Stokeslets, as a method for fluid-structure interaction at zero Reynolds 
number, has been successfully applied to model various biological phenomena. In this application, 
we presented an adaptation of the SBF geometric model developed in (20) to the simulation of open 
planar curves in the method of Regularized Stokeslets, and also introduced a parametric RBF model 
for comparison against the SBF model. 

We conclude the following: 

• SBF and RBF interpolants both show almost identical accuracy when used in the modeling 
of twice-differentiable open shapes. Given the superior accuracy of the SBF interpolant for 
modeling closed periodic shapes, this makes the SBF a promising choice for modeling a 
variety of shapes (closed and open) important to biological applications. 

• The SBF, RBF and Lagrange-Chebyshev interpolants all offer excellent accuracy and 
convergence in the modeling of twice-differentiable open shapes; indeed, as expected, for 
small numbers of points, these interpolants are superior to the Fourier interpolants that are 
sometimes used in the MRS (c/, (20)). 

• The SBF and RBF interpolants are slightly superior to the Lagrange-Chebyshev interpolant 
in terms of accuracy for a given number of data sites. Since SBFs and RBFs likely generalize 
more readily to the modeling of parametric sheets and curves embedded in K 3 , this establishes 
SBFs and RBFs as the more promising candidates for geometric modeling. 

• The modeling of twice-differentiable open curves requires the use of clustered points. In 
the case of Lagrange interpolants, the obvious choice is the set of Chebyshev nodes if the 
clustering is required at the ends of the interval. However, for RBF and SBF interpolants, it 
is clear that the mapped Chebyshev or KTE nodes are superior (for a given shape parameter). 
Indeed, since the KTE nodes used in our work are almost uniformly-spaced, it is reasonable 
to use them as general-purpose node sets for geometric modeling of both closed and open 
curves. 

• The dual representation afforded by global interpolants (data sites and sample sites), when 
used correctly, can offer greater accuracy within the method of Regularized Stokeslets for a 
similar or slightly lower computational cost. This result is true for all interpolants used in this 
work, and for Fourier interpolants as well (which are currently not used in this fashion within 
the MRS). 

When the SBF representation was used within the IB method, a key feature was the use of both the 
data sites and the sample sites (which are analogous to IB points). The number of data sites in that 
work was determined by the statements on geometric accuracy made in (20) , while the number of 
sample sites is chosen such that the spacing (measured in the embedding space) between any two 
sample sites is never more than |, where h is the spacing of the background Eulerian grid. Since the 
method of regularized Stokeslets is grid-free (unlike the IB method), the number of sample sites N s 
is purely determined by the order of the quadrature rule used to compute the velocity (and pressure). 
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The RBF representations (and indeed, other global representations like Fourier-based ones) will 
likely allow the use of higher-order quadrature rules without adversely affecting the accuracy in the 
computation of Lagrangian forces (since the error in the Lagrangian forces converges very rapidly 
with increasing No). This would also allow the use of higher-order time-stepping. We plan to explore 
this possibility in future work. 

We also note that error in the MRS will be due to error in the calculation of Lagrangian forces, 
eiTor from discretizing the structure, or a quadrature rule for Equation 0- as well as error from 
the regularization parameter |2j. In order to have accuracy of force evaluations based on finite 
differences, greater accuracy is achieved provided N s is reasonably large. Thus, the accuracy of the 
pressure and velocity computations depend on the number of IB points N s . As shown in the case of a 
closed curve, we get reduced error for reduced number of points when using parametric interpolants 
to calculate forces. It is also apparent from reviewing the RBF and MRS literature that many of the 
popular regularized blob functions are themselves normalized RBFs. It may be possible to exploit 
this connection to either enable fast evaluation of the forces, velocities and pressures, especially 
if the RBF used for the Lagrangian representation is the same as the RBF that serves as the blob 
function. We hope to explore this connection more carefully in future work. 

Finally, we remark that a non-trivial extension of this current work would be to apply the SBF 
and RBF geometric models to the modeling of centerlines of elastic rods within the method of 
Regularized Stokeslets; since such modeling often requires the application of boundary conditions, 
the SBF and RBF geometric models need to be extended to accomplish this in a rigorous fashion. 
We also plan to explore this in future work. 
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